Dynamical stochastic simulation of complex electrical behavior in neuromorphic networks of metallic nanojunctions

Nanostructured Au films fabricated by the assembling of nanoparticles produced in the gas phase have shown properties suitable for neuromorphic computing applications: they are characterized by a non-linear and non-local electrical behavior, featuring switches of the electric resistance whose activation is typically triggered by an applied voltage over a certain threshold. These systems can be considered as complex networks of metallic nanojunctions where thermal effects at the nanoscale cause the continuous rearrangement of regions with low and high electrical resistance. In order to gain a deeper understanding of the electrical properties of this nano granular system, we developed a model based on a large three dimensional regular resistor network with non-linear conduction mechanisms and stochastic updates of conductances. Remarkably, by increasing enough the number of nodes in the network, the features experimentally observed in the electrical conduction properties of nanostructured gold films are qualitatively reproduced in the dynamical behavior of the system. In the activated non-linear conduction regime, our model reproduces also the growing trend, as a function of the subsystem size, of quantities like Mutual and Integrated Information, which have been extracted from the experimental resistance series data via an information theoretic analysis. This indicates that nanostructured Au films (and our model) possess a certain degree of activated interconnection among different areas which, in principle, could be exploited for neuromorphic computing applications.

ij , ∀ij Require: rnd is a uniformly distributed random number, ∈ [0, 1) for ij ∈ network do loop over all netwrok's edges if σ ij = σ α then if (|∆V ij | > ∆V th ) and (rnd < P nl ) then σ ij ← σ β try to apply nonlinear update move else try to apply thermal dissipation-based update move (a link with σ ij = σα can only be upgraded) ij ≥ W th 0 up ) and (rnd < P up ) then ij ≥ W th 1 up + W th 0 up ) and (rnd < P up ) then ij > W th 2 up + W th 1 up + W th 0 up ) and (rnd < P up ) then σ α ← σ δ else if σ ij = σ β then try to apply nonlinear update move if (∆V th > |∆V ij |) and (rnd < P nl ) then σ ij ← σ α else try to apply thermal dissipation-based update move ij ≥ W th 1 up ) and (rnd < P down ) then try to apply thermal dissipation-based update move ij ≥ W th 2 down ) and (rnd < P down ) then σ ij ← σ γ try to apply thermal dissipation-based update move (a link with σ ij = σ δ can only be downgraded) Here follows a list of optimal parameters for the simulation input and for the MC moves.
• thermal dissipation based moves: P up = P down = 0.0015 (probability to accept any move that changes σ ij upwards or downwards, regardless of the specific value of σ ij ): An example will immediately clarify the link update mechanism, according to its thermal down , the link conductance will be decreased down to σ β = 10 −3 1/Ω with probability 0.0015, and so on. A completely identical mechanism holds for the absorbed power coming from each link neighbors. Obviously, links whose initial conductance is σ ij = σ α = 10 −11 1/Ω can not be further downgraded.
Clearly, these particular choices of the threshold and probability values are not unique, but this set allows to approximately retrieve the essential experimental features in the simulated data, as largely discussed in the paper. Slight variations of these thresholds were found to generate analogous results. Figure S1 represents I(∆V ) cycles performed with networks of different sizes, aiming to highlight the minimal complexity required for the qualitative reproduction of the experimental data. In these simulations, ∆V ranges from -8 V to 8 V, with 8000 MC steps for each applied voltage, and each network size has been investigated with 10 independent runs. If the network features a single layer, made by N x × N y = 32 × 15 nodes, the system presents a nonphysical behavior with wild current oscillations (see top panel). The errorbars are almost anywhere extremely large, and a similar tendency is retrieved even when a second layer is stacked on the top of the first one (panel b). This can be ascribed to the small amount of links in the network which, in turn, induce huge fluctuations in the availability of effective paths for the current to flow through the system. The addition of a third layer induces the onset of two distinct regimes, at low and high voltages, with a partial reduction of the errorbars (bottom panel). Fig. S1c begins to feature some similarities with the analogue experimental picture. We decided to use 3 layers and a much larger network (N x × N y = 42 × 27), with the aim to further reduce the errorbars onto the I(∆V ) curves and to provide more alternative paths for the input current to flow towards the output node. Nonetheless, the complexity of our system clearly remains much smaller than the experimental sample one, and this is -for instance -the origin of the current saturation at high voltages which is present in our SRN model, even in the largest simulated one.

SUPPLEMENTARY NOTE: SHORTEST PATHS ANALYSIS
As mentioned in the main text, the structure and the topology of the network were studied via the investigation of the shortest path connecting the source and the sink of the system, thus exploring the possible paths in which the current could flow. This has been achieved thanks to the analysis we carried out with NetworkX tools [44]. We anticipate here that the application of a high voltage is found to dramatically reduce the number of alternative options to reach the output node.
What does it mean 'shortest'? Distances between pairs of nodes within a graph can be indeed measured, once that a metrics has been chosen. We consider here the I matrix (N n × N n , i.e. 3404 × 3404) as a weighted graph for each MC step. Here we sketch the procedure, based on a Python custom code which integrates some NetworkX functions: • load, for each time frame t saved during the simulation, the corresponding matrices I| t and A| t into 2D numpy arrays.
• loop over I| t elements: each entry I ij of I| t corresponds to a weight The current values cover a wide range, typically between 10 −11 A and 10 −2 A. The rationale behind the construction of this weight matrix W is that those links where the conductance is non-vanishing should be used by the current, while connecting node 0 and node N n − 1. Conversely, highly (infinitely) resistive edges are essentially never crossed by the electrical current; assigning them an infinite weight means excluding these edges (because the nodes they link are substantially disconnected) from the search of the SP connecting the source and sink nodes.
• once built the weight matrix W, given two nodes i and j, it is possible to calculate the shortest path between them, where the notion of distance is given by the weight matrix W.
NetworkX includes some methods for calculating the shortest path between a pair of nodes, in particular those by Dĳkstra [10] and Bellman-Ford [11,12].
Therefore, two quantities of interest can be defined: • L SP , i.e. the length of the shortest path, measured as an inverse current in 1/A. We underline that, with this metrics, the shortest paths correspond to the most conductive source-drain pathways.
• N SP , the number of equivalent shortest paths In Figs. S2 and S3 we show, for a typical simulation, L SP · ∆V (which has the units of a resistance) and of N SP as a function of the MC step t M C . On the one hand, it is evident that the length of the shortest path is strongly dependent on the applied voltage: Figure S2 shows the evolution of L SP · ∆V for a simulation with ∆V tot = 1 V (grey), before the application of the high voltage phase (green), which (at equilibrium, around 10000 MC steps) is very similar to the equilibrium L SP reached by the system at the same voltage after the conditioning step (blue).

SUPPLEMENTARY NOTE: SPATIAL COARSE-GRAINING
In Fig. S4 it is shown the spatial coarse-graining method adopted for the network analysis via Information Theory tools. Black lines are the boundaries of a parallelepiped whose links are sent onto the coarse-grained region with index 7, as an example. An analogous procedure sends the other 6 subregions of the original system (left), onto the corresponding coarse-grained indexes (right). Note that also links connecting source/sink at the nods of the first layer participate into the spatial coarse-graining. Links of the original systems are colored according to their conductance, apart from the σ α ones which are left white.
FIG. S4. Left: a typical snapshot of the resistor network, with link conductances colored with the same color scale used in Fig. 4. Realized with NetworkX [9]. Right: scheme of the spatially coarse-grained system obtained by dividing the network in 7 sub-regions. The sub-region of the original system which is mapped onto the coarse-grained group of index 7 is bounded by thick black lines.